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Abstract 

This paper reviews the Serret-Andoyer (SA) canonical formahsm 
in rigid-body dynamics, and presents some new results. As is well 
known, the problem of unsupported and unperturbed rigid rotator 
can be reduced. The availability of this reduction is offered by the 
underlying symmetry, that stems from conservation of the angular 
momentum and rotational kinetic energy. When a perturbation is 
turned on, these quantities are no longer preserved. Nonetheless, the 
language of reduced description remains extremely instrumental even 
in the perturbed case. We describe the canonical reduction performed 
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by the Serret-Andoyer (SA) method, and discuss its apphcations to 
attitude dynamics and to the theory of planetary rotation. Specifi- 
cally, we consider the case of angular-velocity-dependent torques, and 
discuss the variation-of-parameters-inherent antinomy between canon- 
icity and osculation. Finally, we address the transformation of the 
Andoyer variables into action-angle ones, using the method of Sadov. 

1 INTRODUCTION 

One of the classical problems of mechanics is that of a free motion of a rigid 
body, usually referred to as the Euler-Poinsot problem. The formulation 
and solution of this problem are usually performed in two steps. First, the 
dynamics of the rotation are represented by differential equations for the 
components of the body angular velocity. As is well known, these equa- 
tions admit a closed-form solution in terms of Jacobi's elliptic functions pLj. 
Second, the kinematic equations are utilised to transform the body angu- 
lar velocity into a spatial inertial frame. While this classical formulation is 
widespread among engineers, astronomers exhibit a marked preference for 
formulation that takes advantage of the internal symmetries of the Euler- 
Poinsot setting [2j. The existence of internal symmetries indicates that the 
unperturbed Euler-Poinsot problem can be reduced to a smaller number of 
variables, whereafter the disturbed setting can be treated as a perturbation 
expressed through those new variables (the disturbance may be called into 
being by various reasons - physical torques; inertial torques emerging in non- 
inertial frames of reference; non-rigidity of the rotator). 

There exist reasons for performing this reduction in the Hamiltonian form 
and making the resulting reduced variables canonical. The first reason is 
that a perturbed Hamiltonian system can be analytically solved in any order 
over the parameter entering the perturbationj^ A celebrated example of a 
successful application of the Hori-Deprit method to a geophysical problem 
is given by the theory of rigid Earth rotation offered by Kinoshita ^ and 
further developed by joint efforts of Escapa, Getino, & Ferrandiz [9l [TO] . 
Getino & Ferrandiz [11], and Kinoshita & Souchay [I2]. Other applications 

^ While in the past such analytical solutions were built by means of the von Zeipel 
method, more efficient is the procedure independently offered by Hori [3] and Deprit [1]. 
An explanation of the method can be found in [5]. A concise introduction into the subject 
can be found in Chapter 8 of the second volume of [6] and in Chapter 5 of 0- 
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to planetary and liunar rotation may be found in [T7], [IB], and 

The second advantage of the canonical description originates from the 
convenience of numerical implementation: symplectic schemes are well-known 
for their good stability and precision. (See, for example, Yoshida [I3].) This 
is why the Hamiltonian methods are especially beneficial at long time scales, 
a feature important in astronomy. (Laskar & Robutel [13]; Touma & Wisdom 

mm) 

Several slightly different sets of canonical variables are used for modeling 
rigid-body dynamics and kinematics. Most popular is the set suggested in 
1923 by Andoyer [20]. This set is not completely reduced: while three of 
its elements are constants (in the unperturbed free-spin case), three others 
are permitted to evolve. Andoyer arrived at his variables through a pretty 
sophisticated procedure. Canonicity of Andoyer's change of variables is not 
immediately apparent (though can be proven by a direct construction of the 
corresponding generating function). Much later, the study by Andoyer was 
amended by Deprit [21j who demonstrated that the canonicity of Andoyer's 
transformation may be proven by using differential forms and without resort- 
ing to explicitly finding a generating function. 

A full reduction of the Euler-Poinsot problem, one that has the Andoyer 
variables as its starting point, was recently offered by Deprit & Elipe [2]. 
It would be interesting to notice that, historically, the pioneer canonical 
treatment of the problem, too, was formulated in a completely reduced way, 
i.e., in terms of canonically-conjugate constants of motion. These constants, 
presently known as Serret elements, are named after the 19th century French 
mathematician Joseph Alfred Serret who discovered these variables by solving 
the Jacobi equation written in terms of the Eulerian coordinates [23]. Serret 's 
treatment was later simplified by Radeau Pl] and Tisserand [2S]. However, 
for the first time the Serret elements appeared in an earlier publication by 
Richelot [26]. 

Probably the most distinguishing feature of the canonical approach is that 
it permits reduction of the torque-free rotational dynamics to one-and-a-half 
degrees of freedom. In essence, such a formulation reduces the dynamics 
by capturing the underlying symmetry of the free rigid body problem - the 
symmetry taking its origin from the conservation of energy and angular mo- 
mentum. As a result, the entire dynamics can be expressed by differential 
equations for two of the Eulerian angles and one of the conjugate momenta 
- equations that are readily integrable by quadrature. The corresponding 
one-and-a-half-degrees-of- freedom Hamiltonian yields a phase portrait, that 
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is similar to that of the simple pendulum and contains a separatrix confining 
the librational motions [2T] . 

In the canonical formulation of attitude-mechanics problems, incorpo- 
ration of perturbation and/or control torques into the picture is a subtle 
operation. In particular, when the perturbing torque is angular-velocity- 
dependent, the canonicity demand contradicts the osculation condition. In 
other words, the expression for the angular velocity through the canonical 
Serret or Andoyer variables acquires a correction called the convective term. 

To better understand the latter observation, we shall use an orbital- 
dynamics analogy. This shall be convenient, because both orbital and ro- 
tational dynamics employ the method of variation-of-parameters to model 
perturbing inputs. In both cases, the coordinates (Cartesian, in the orbital 
case, or Eulerian in the rotational case) are expressed, in a non-perturbed 
setting, via the time and six adjustable constants called elements (orbital 
elements or rotational elements, accordingly). If, under disturbance, we use 
these expressions as an ansatz and endow the "constants" with time de- 
pendence, then the perturbed velocity (orbital or angular) will consist of a 
partial derivative with respect to the time, plus the convective term, one that 
includes the time derivatives of the variable "constants." Out of sheer conve- 
nience, the so-called Lagrange constraint is often imposed. This constraint 
nullifies the convective term and, thereby, guarantees that the functional 
dependence of the velocity upon the time and "constants" stays, under per- 
turbations, the same as it used to be in the undisturbed setting. The variable 
"constants" obeying this condition are called osculating elements. Otherwise, 
they are simply called orbital elements (in orbital mechanics) or rotational 
elements (in attitude mechanics). 

When the dynamical equations, written in terms of the "constants," are 
demanded to be canonical, these "constants" are the Delaunay elements in 
the orbital case, or the initial values of the Andoyer elements, in the spin 
case. These two sets of elements share a feature not readily apparent: in 
certain cases, the standard equations render these elements non-osculating. 

In attitude dynamics, the Andoyer variables come out non-osculating 
when the perturbation depends upon the angular velocity. For example, since 
a transition to a non-inertial frame is an angular-velocity-dependent pertur- 
bation, then amendment of the dynamical equations by only adding extra 
terms to the Hamiltonian makes these equations render non-osculating An- 
doyer variables. To make them osculating, extra terms must be added in the 
dynamical equations (and then these equations will no longer be symplectic) . 
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Calculations in terms of non-osculating variables are mathematically valid, 
but their physical interpretation is not always easy. 

The purpose of this paper will be threefold. We shall first aim to provide 
the reader with a coherent review of the Serret-Andoyer (SA) formalism; sec- 
ond, we shall dwell on the canonical perturbation theory in the SA context; 
and third, we shall consider introduction and physical interpretation of the 
Andoyer variables in precessing reference frames, with applications to plan- 
etary rotation. A subsequent publication by Bloch, Gurfil, & Lum (2007) 
[22j will deal with the Andoyer-variables-geometry based on the theory of 
Hamiltonian systems on Lie groups, and will fill the existing gap in the con- 
trol literature by using the SA modeling of rigid-body dynamics in order to 
derive nonlinear asymptotically stabilising controllers. 

2 THE EULERIAN VARIABLES 

Through this entire section the compact notations = sin(x) , = cos(x) 
will be used. The angular-velocity and angular-momentum vectors, as seen 
in the principal axes of the body, will be denoted by low-case letters: a; and 
g . The same two vectors in an inertial coordinate system will be denoted 
with capital letters: fl and G . 

2.1 Basic formulae 

Consider a rotation of a rigid body about its center of mass, O. The body 
frame, B, is a Cartesian dextral frame that is centered at the point O and is 
defined by the unit vectors hi, b2, constituting the fundamental plane, and 
by = bi X b2. The attitude of B will be studied relative to an inertial 
Cartesian dextral frame, /, defined by the unit vectors Si, §2 lying on the 
fundamental plane, and by S3 = Si x §2. These frames are depicted in Fig. [2j 
A transformation from I to B may be implemented by three consecutive 
rotations making a 3 — 1 — 3 sequence. Let the line of nodes (LON) OQ2 be 
the intersection of the body fundamental plane and the inertial fundamental 
plane, as shown in Fig. [2l Let 1 be a unit vector pointing along 0Q2- The 
rotation sequence can now be defined as follows^l 

^ Be mindful that in the physics and engineering hterature the Euler angles are tra- 
ditionally denoted with ((j) , 6 , ip) . In the literature on the Earth rotation, the inverse 
convention, ( -0 , 6* , ) , is in use. P [9l [TO l [T H [T2] 
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• -R(0, S3), a rotation about §3 by < < 27r, mapping Si onto 1 ; 

• R{9^ i), a rotation about 1 by < < vr, mapping §3 onto b3 ; 

• R{il),h'i), a rotation about bs by < < 27r, mapping 1 onto bi . 

The composite rotation, R G 50(3), transforming any inertial vector into 
the body frame, is given by 



R{ct),e,i,) = R{ijM)R{e;\)R{ct),s^) 
Evaluation of this product gives 



(1) 



i?(0,^,^) 



c^C0 - s^ceS0 c^s,^ + s^cgc,^ s^s^ 
— s^C(^ — c^cqs^ —s^s^ + c^cec^ c^se 

SeS^ -SgC^ Cg 



(2) 



To write the kinematic equations, we recall that the body angular-velocity 
vector, uj = [cui, a;2, Ws]"^ satisfies 



where the hat map 



-RR^ (3) 
so (3) is the usual Lie algebra isomorphism and 



U!2 

LU^ — Co"! 

—UJ2 uji 



(4) 



Insertion of Eq. ([2]) into ([3]) yields the well-known expressions for the com- 
ponents of the vector of the body angular velocity uj : 



uji = (/)Ses^ + 9c^ , 

UJ2 = 0C^S6I - 6s^ , 

U3 = i' + (j)Ce , 



(5) 
(6) 
(7) 



while the action of the rotation matrix ([2]) upon the body angular velocity 
gives the inertial-frame-related (sometimes called spatial) angular velocity: 



R^ UJ 



(8) 
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with the following components: 

Qi = 9 + tp se , (9) 
fl2 = 6 + ^ se , (10) 
= ^ + ^ce . (11) 

In the body frame, attitude dynamics are usually formulated by means of 
the Euler-Poinsot equations. In a free-spin case, these equations look as 

Id; + a; X lo; = , (12) 

I being the inertia tensor. To arrive to these equations, one has to start out 
with the Lagrangian, 

/:(0,^,V^,0,^,^) = ia;-Ia; , (13) 

to substitute therein expressions (12121 - [2T4l) : to write down the appropriate 
Euler-Lagrange variational equations for the Euler angles; and then to use 
the expressions (12121 - 12 141) again, in order to rewrite these Euler-Lagrange 
equations in terms of uj . This will result in (1121) . The same sequence of 
operations carried out on a perturbed Lagrangian £ + A£ will produce the 
forced Euler-Poinsot equations: 

Icj + a; X Icj = u , (14) 

u G M'^ being the body-frame-related torque that can be expressed through 
derivatives of A£ . The Euler-Poinsot description of the motion is essentially 
Lagrangian. Alternatively, the equations of attitude dynamics can be cast 
into a Hamiltonian shape. 

A most trivial but important observation can be made simply from look- 
ing at f|T4l) . If one rewrites f|T^ not in terms of velocities but in terms of 
the Euler angles, three differential equations of the second order will emerge. 
Their solution will depend on the time and six adjustable constants. Hence, 
no matter which description one employs - Lagrangian or Hamiltonian - the 
number of emerging integration constants will always be six. 

2.2 Hamiltonian description. The free-spin case. 

From here forth we shall assume that the body axes coincide with the prin- 
cipal axes of inertia, so we can write 

I = diag(Ji,/2,/3). (15) 
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Having substituted expressions (12121 - 12141) and (fT^ into (fT^ , one can eas- 
ily write the generalised momenta Pn = $ , 6 , \E' conjugate to the 
configuration variables = cj) , 9 , ip [21 [22] : 

dC ■ ■ ■ ■ ■ ■ 

$ = — ^ = /iS0S^(<^s^S0 + 6'c^) + /2Sec^(</)c^se - 6's^) + /3Ce(</)c6' + ^Z)) , (16) 

0(p 

dC ■ ■ ■ ■ 

= — X = hc^{4>S0S^ + ec^) - /2S^((/)C0Se - 6s^) , (17) 

o9 
dC 

VP = =/3(0c, + ^) . (18) 

dip 

The inverse relations will assume the form of 

• _ /iC^(^'C6lC^ - $C^ + QsgS^) + /2S^(^C0S^ - $S^ - Oc^S^) 



■ /2C^($S^ + GsgC^ - ^S^Cg) + /iS^(^C^C9 - $c^ + Osgs^) 

" = Ja^. ■ ""^ 

^ _ /i/3C^('l>cec^ - \&c^Cg - 8s^secg) + /3/2S^('l>ces^ - \&c^s^ + Gc^sgCg) 

-^1-^2-^380 

+ Y (21) 
Substitution of f[T^ - I^Tl) into the Legendre-transformation formula 

n = ^^ + ee + -^^-c (22) 

will then readily give us the free-spin Hamiltonian: 



+ ,l_iUi^,es,c,. (23) 



In the above Hamiltonian the coordinate (j) (the angle of rotation about the 
inertial axis S3 ) is cyclic (ignorable), so that the appropriate momentum 
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$ is an integral of motion. This symmetry implies a possibility of reduction 
of the free-spin problem to only two degrees of freedom. Based on this 
observation, Serret |23j raised the following question: Is there a canonical 
transformation capable of reducing the number of degrees of freedom even 
further? We shall discuss this issue in the following section. To get this 
program accomplished, we shall need to know the relationship between the 
body rotational angular momentum and the conjugate momenta. Let g = 
Oi^i and G = ^ GiSi be the angular momentum in the body frame and 
in inertial frame, accordingly. By plugging (12121 - [2T^ into 

g = Icv (24) 

and by subsequent insertion of (fTOl - [21]) therein, one easily arrives at 

$S^ + QsgC^ - "^S^Cg 

91 — , l^oj 

$C^ - Qs^Sg - ^C^Ce 

92 = — , (26) 

sg 

93 = . (27) 

Notice the symplectic structure defined by the components of the angular 
momentum. It is a matter of computing partial derivatives to check that the 
Poisson brackets are 

{91,92) = -93, {92,93) = -9i, {93,91) = -92 (28) 
By the same token, insertion of expressions ([9] - [TTj) into 

G = m , (29) 
with the subsequent use of (fT9]- [2T!). entails: 

^ _ ^S0 + QsgC^ - $S0Cg ^g^^ 

^CgC^ - "^C^ + Qs^Sg 
Gr2 — , [oi) 

G3 = $ , (32) 
yielding the symplectic structure 
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(Gi,G2) — Ca, {G2,G3) — Gi, (Gs, Gi) — (j2- (33) 

from the symplectic structure (!28l) and the expression (150|) of the Hamilto- 
nian, there results 



1 
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gi = {g,,n) = -\^---jg293, (34) 
92 = {92,'H) = - ( ^ - ^] gsgi, (35) 



^1 = {92,T-C) = - yj-- —j gig2. (36) 

This system is integrable because it admits two integrals, the energy ([T3|) 
and the norm of the angular momentum |G| = G. With this integral, we 
may regard the phase space of flMl) - fl36|) as a foliation of invariants manifolds 

S'iG) = {igi,g2,93M + 9l + 9l = G'}. (37) 
The trajectories will be the level contours of the energy ellipsoid, 



T..„ = i-.G = i(/,..? + /2c.,^ + /3..3^) = i(| + | + |) 



(38) 



on the sphere (l37|l . as can be seen in Figure [H 

Below we shall need also the expressions for the canonical momenta via 
the components of the angular momentum G. As agreed above, the Euler 
angles 0, 6, ip determine the orientation of the body relative to some inertial 
reference frame. Let now the angles (po, J, I define the orientation of the body 
relative to the invariable plane (one orthogonal to the angular-momentum 
vector G), as in Fig. [21 and let h, I, g — (po be the Euler angles defining the 
orientation of the invariable plane relative to the reference one|§ Evidently, 

gi = G sj si , g2 = G sj ci , g^ = G cj . (39) 

It is now straightforward from (12^ - 1271) and fl5^ that 

^ = 9ises^ + g2 sgc^ + g^ce = G ( ce cj + sg sj C(^_/) ) = G cj , (40) 



It would be natural to denote the orientation of the body relative to the invariable 
plane with 0o , , V'o , but we prefer to follow the already established notations. 
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Q = giC^ - 923^ = G Sj S(i_^) 



(41) 



m = = Gcj 



(42) 



3 THE SERRET-ANDOYER 
TRANSFORMATION 

3.1 Richelot (1850), Serret (1866), Radau (1869), 
Tisserand (1889) 



The method presently referred to as the Hamihon-Jacobi one is based on the 
Jacobi equation derived circa 1840. Though Jacobi's book [33] was pubhshed 
only in 1866, the equation became known to the scientific community already 
in 1842 when Jacobi included it into his lecture course. In 1850 Richelot 
suggested six constants of motion grouped into three canonically-conjugate 
pairs. These constants became the attitude-dynamics analogues of the Delau- 
nay variables emerging in the theory of orbits. Later Serret [23] wrote down 
the explicit form for the generating function responsible for the canonical 
transformation from {ip, 6, 0, 0, $) to Richelot's constants. His treatment 
was further polished by Radau [2^ and explained in detail by Tisserand [25] 
The canonical transformation, undertaken by Serret, 



( gi = (p, q2 = 0, qs = ij , pi = $ , p2 = , P3 = ^ ; ^(g, p) ) 



is based on the fact already known to the mathematicians of the second 
part of the XIX*^ century: since both descriptions, ( q, p, 7Y(q, p) ) and 
(Q, P, 7i*(Q,P)) , are postulated to satisfy the Hamiltonian equations, 
then the infinitesimally small quantities 



{ Ql , Q2 , Q3 , -Pi ; -P5 



'2, A; n*iQ,P)) 



(43) 



p dq^ — H dt 



(44) 



and 



dC = Cf dP + n* dt 



(45) 
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are perfect differentials, and so is their sum 



dS = dC + dC= p^rfq + Q'^ dP - (n - W) dt . (46) 

If we start with a system described with ( q, p, ?^(q, p) ) , it is worth looking 
for such a re-parameterisation (Q, P, 7i*(Q,P)) that the new Hamilto- 
nian H* is constant in time, because this will entail simplification of the 
canonical equations for Q and P . Especially convenient would be to find a 
transformation that would nullify the new Hamiltonian Ti* , for in this case 
the new canonical equations would render the variables ( Q, P ) constant. 
One way of seeking such transformations is to consider 5* as a function of 
only q , P , and t . Under this assertion, the above equation will result in 

f)Q f)Q f)Q 

dt + —dq+ —dP = p^dq + Q^dP - {n- n*) dt (47) 
at oq oP 

whence 

^ - dP ' m - ^ ■ ^^^^ 

The function ^(q , P , t) can be then found by solving the Jacobi equation 

dS \ dS ^,,fdS „ \ 
niq,—,t)+ — = n*i—,P,t) . (49) 



dq' J dt \dP 
In the free-spin case, the Jacobi equation becomes 

1 / sin^ cos^ ij\ f dS dS ^\ ^ 1 / cos^ sin^ 'tb\ f dS\'^ 



2smH\ h h J \d(p di) J 2 \ h h J \d9 J 



1 fdSV fl l\dS fdS dS \ sin^ cos^ OS ^,,fdS ^ \ 

At this point, Serret [23], Radau [21], and Tisserand [25J chose to put Ti* 
equal to zero, thereby predetermining the new variables ( Q , P ) to come 
out constants. (See, for example, equation (21) on page 382 in [25].) Thence 
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the Jacobi equation (HUl) became equivalent tc0 

dS\ dS dS dq 

By taking into account that the initial Hamiltonian Ti depends explicitly 
neither on the time nor on the angle (f) , Serret and his successors granted 
themselves an opportunity to seek the generating function in the simplified 
form of 

the constant Ai being equal to the negative value of the Hamiltonian Ti , 
i.e., to the negative rotational kinetic energy: 

M = - Tkin ■ (53) 

The second constant, A2 , as well as the derivatives dS/dO and dS/dip 
can be calculated via the first formula (HHI): 



SO we get: 

5 = - t Tfci„ + $ + [ Ode + I ^ d^j + C , (55) 



$ , G , and ^ being given by formulae (H0]- H2|) . and $ being a constant of 
motion because Ti is 0-independent. Plugging of (jlO-llZ]) into ( 155|) yields: 

S = - t Tkin + $ + G j COS J dip + G j sinJ sm{l - tp) de + C 
^ The generating function should be written down exactly as ([5T|) and not as 

because the new momenta P arc not playing the role of independent variables but will 
emerge as constants in the solution. 
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= — t Tkin + G (f) COS I + G COS J dl + 



G j [- cosJ d{l - tp) + sin J sin {I - ^j) d9] + G , (56) 

The latter expression can be simplified through the equality (derived in the 
Appendix A.l) 

d g — d {(j) — h) cos I = — cos J d{l — ip) + sin J sin (/ — ip) dO (57) 

the angles (0 — h) and g being shown on Fig. [2l and Fig. [31 This will 
entail: 

S = - tTkin + Gh cos J + G g + G j cos J dl + G , (58) 

where the constant G may be chosen, for example, as zero. At this point, 
the authors of the XIX*'^ century made use of two observations. First, in the 
unperturbed case, one can introduce the dimensionless time 

u = n {t - to) , (59) 

n = ^(G^/f^ — 2Tkm) {l2^ ~ ^3^) being the mean angular velocity 
(or, playing an astronomical metaphor, the "mean motion"). At each in- 
stance of time, the value of so-defined u is unambiguously determined by 
the instantaneous attitude of the rotator, via the solution of the equations of 
motion (this solution is expressed through the elliptic functions of u). More 
specifically, m is a function of T^in , G , G cos J , 9 , and ip - the appropri- 
ate expression is given in Appendix A. 2 below. Simply from this definition 
it trivially follows that in the unperturbed case the difference u/n — t is an 
integral of motion, because it is identically nil. Under perturbation, though, 
the afore mentioned function uiT^in , G , G cos / , 9 , ip) , divided by n , will 
no longer furnish the actual values of the physical time t . (It would, had we 
substituted in fl59|) the product with an integral.) Hence, under perturbation, 
the integral of motion u/n — t will acquire time dependence. 

The second observation can be done if we consider a coordinate system 
defined by the angular-momentum vector and a plane perpendicular thereto. 
This plane (in astronomy, called invariable) will be chosen as in Fig. [2] and 
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Fig. [31 so that the centre of mass of the body he in this plane. The Euler 
angles defining the attitude of a body relative to the invariable plane are 
g , J , / , as on Fig^ [2l and Fig. [31 In the case of free spin, they obey the 
differential relation! 



J5 



dg + cos J dl = — dt (60) 
whose integration results in 

9 - - + / cos J dl = go . (61) 

G n J 

Trivially, the expression standing in the left-hand side is an integral of motion 
in the unperturbed case, because in this case it is equal to go ■ However, 
under disturbance the left-hand side will no longer be equal to go (and thus 
will no longer remain an integral of motion), because neither fISI?]) nor fpUl) 
will hold. 

The above two observations motivate the following equivalent transfor- 
mation of f[58l): 



S = (2 - Tkin + Gh cos/ + Ggo . (62) 

Now it might be tempting to state that, if we choose (-Pi, P2 , -P3 ) = 
( G , G cos / , Tkin ) , then expression (!62l) would entail 

and 

dS 



d{G cos/) 



Q2 (64) 



^ Tisserand (PS], p. 386) and his contemporaries used to prove (|M)) by a somewhat 
tedious method based on formulae connecting the angles J , I with the components of 
the body angular velocity u: and with the absolute value of the angular momentum: 

IiLUi — G sin J sin I , I2W2 = G sin J cos I , I3UJ3 — G cos J . 

We provide this proof in the Appendix A. 3 below. 

The Andoyer variables, which we shall introduce below, make (|60|) self-evident: (|60|) 
follows from equations ([921), jOSl) and dM]). 
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where 



Qi = 9o (65) 

and 

Q2 = h (66) 

are integrals of motion. These two formulae, if correct, would implement our 
plan set out via the second equation 

In fact, equations fl62]-[63l) do not immediately follow from fl6T|) . because 
the variations of G , G cos / , Tkm are not independent from the variations 
of the angles involved. These variations are subject to some constraints. 
A careful calculation, presented in the Appendix below, should begin with 
writing down an expression for 6S , with the said constraints taken into 
account. This will give us: 

55 = (^^ - 6Tkin + h6{G cos/) + go 6G (67) 



whence we deduce that fl63] - 1661) indeed are correct. We also see that the 
negative of the initial instant of time (which, too, is a trivial integral of 
motion), turns out to play the role of the third new coordinate: 



dS 



dTi 



Q3 (68) 



kin 



where 



u 

= - - t = - to . (69) 
n 

All in all, we have carried out a canonical transformation from the Euler 
angles and their conjugate momenta to the Serret variables Qi = go , 
Q2 = h , Q3 = — to and their conjugate momenta Pi = G , P2 = 
G cos/, P3 = Tkin ■ In the modern literature, the constant P2 is denoted 
by H . Constants h , G , H enter the Andoyer set of variables, discussed 
in the next subsection. 

3.2 Andoyer (1923), and Deprit & Elipe (1993) 



While Serret [2^ had come up with a full reduction of the problem, Andoyer 
suggested a partial reduction, based on a general method of obtaining 
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a set of canonical variables {n coordinates and their n conjugate moments) 
from a set of 2n variables that are not necessarily canonical. This general 
method is discussed in the Appendix, §A.5. 

To understand the essence of that partial reduction, we shall start with a 
transformation from the inertial to the body frame via a coordinate system 
associated with the invariable plane. 

Referring to Fig. [21 let OQi and OQ3 denote the LON's obtained from 
the intersection of the invariable plane with the inertial plane and with a 
plane fixed within the body, respectively. Let i be a unit vector along the 
direction of OQi, and j be a unit vector along OQ3. Define a3 — 1 — 3 — 1 — 3 
rotation sequence as follows: 

• R{h, S3), a rotation about S3 by < /i < 27r, mapping Si onto i. 

• R{I,i), a rotation about i by < / < vr, mapping §3 onto the angular 
momentum vector, G. 

• R{g,G/G), a rotation about a unit vector pointing in the direction of 
the angular momentum hj < g < 2tt, mapping i onto j. 

• i?(J, j), a rotation about j by < J < vr, mapping G onto b3. 

• R{l,hs), a rotation about b3 by < / < 2n, mapping j onto bi. 
The composite rotation may be written as 



R{h, I, g, J, I) = R{1, h3)R{J,3)R{g, G/G)R{I, i)R{h, S3) 



(70) 



Evaluation of the product of these five matrices will result in 



R{h,I,g, J,l) = [vi, V2, V3] 



(71) 



where 



{ClCg - SlCjSg)Ch - {ClSg + S,CjCg)C/ " SlSjSjSh 
-S/Cg - ClCjSgCh - (-S/Sg + C/CjCg)C/ " C/SjS/S,, 

SjSgCh + (sjCgC/ + CjS/)s,, 



(72) 



(C/Cg - SlCjSg)Sh + (C/Sg + SiCjCg)C/ - S/SjS/Cft 

-SiCg - CiCjSgSh + (-S/Sg + C/CjCg)c7 " C/SjS/C/j 
SjSgSh + (-SjCgC/ - CjS/)c,, 



(73) 
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(C/Sg + SiCjCg)si + S/SjC/ 

V3 = -siSg + ciCjCgSj + cisjcj (74) 

-SjCgS/ + CjC/ 

A sufficient condition for the transformation fH3|) to be canonical can be 
formulated in terms of perfect differentials. Given that the Hamiltonian lacks 
explicit time dependence, this condition will reaci§ [50] : 

$^0 + Qde + = Ldl + Gdg + Hdh (75) 

Let us ffist evaluate the left-hand side of (!75l) . The differential of R{(j),9,ip) 
is readily found to be [35] 



dR{(f), 9, ip) = s^dcj) + \dd + b3# 



(76) 



Multiplying both sides of (1751) by G, and taking advantage of the identities 
(cf. fl30]- [32]) ). one will arrive at 



G ■ S3 

whence 



G ■ i = G ■ (sic^ + S2S<^) = e , G 



^ , (77) 



^d<p + Qde + ^dtjj = 



G ■ S3 rf0 + G ■ (si C0 + §2 s^) d9 + G ■ h-s dip = G ■ dR . (78) 

We shall now repeat the above procedure for the right-hand side of (!75|) . The 
differential of i?(/i, /, g, J, /) is evaluated similarly to (1761) so as to get 

dR{h , I , g, J , I) = S3 dh + I dl + — dg + j dJ + dl (79) 
Since 

G-S3 = $, G-i = 0, G-G/G' = G', G-j = 0, G ■ 103 = * , (80) 
then multiplication of both sides of fl79l) by G will lead us to 

G-dR= Mh + Gdg + ^dl (81) 

^We shall show shortly that indeed the perfect differentials criterion for canonical trans- 
formations holds in this case with the perfect differential of the generating function being 
identically equal to zero. 
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Together, ([7S]) and ([HI]) give 



(l> d(p + e d9 + dij = G ■ dR = <i> dh + G dg + ^ dl , (82) 

comparison whereof with fl75p immediately yields 

^ = H , ^ = L . (83) 

To conclude, the transition from the Euler coordinates (f) , 6 , to the coor- 
dinates h , g , I becomes a canonical transformation (with a vanishing gen- 
erating function, as promised above), if we choose the momenta H , G , L 
conjugated to h, g , I as H = ^ , G = \G\ , L = "if , correspondingly. 

To get direct relations between the Euler angles and Andoyer angles 
h, g, / , we can compare the entries of R{(l),9,ip) and R{h, I, g, J,l). Al- 
ternatively, we may utilise the spherical laws of sines and cosines written for 
the spherical triangle Q1Q2Q3 [8]- Adopting the former approach, we first 
express / and J in terms of the SA canonical momenta (cf. Fig. [2]), 

ci = H/G , cj = L/G . (84) 

Equating the (3, 3) entries, the quotient of the (1, 3) and (2, 3) entries, and 
the quotient of the (3, 1) and the (3, 2) entries yields, respectively, 

T H c 

= - ^V(G2-L2)(G2-if2) (85) 



. V^^^s,c. + L.c,)+i/.v/G^ 
-VG^ - H\G^gSi + Lcicg) + Hci^G^ - 



tan (A = - L\GsgCn + /7s,c,) + LsfeVG^ - 

VG^ - L2(Gs,s;, - HchCg) - LchVG^ - 

To complete the equations of transformation relating the Eulerian vari- 
ables to the Andoyer variables, we utilise either Eqs. fl90l) - fl92l) or Eqs. fl30|) - 
f l32|) to calculate the magnitude of G and then express the momentum in 
terms of the Andoyer variables. Carrying out this procedure yields 

= Gsjs^,-; = - LH^_i (88) 
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where the expressions for and s^ may be readily obtained from Eq. 

Substituting the identities (1551) and (15^ - [HH]) into ([23]) yields the new, 
single-degree-of-freedom Hamiltonian 

n{g, h, /, G, ^) = m + I) - + ^ 
Alternatively, we could have utilised relations (!39l) and (!84l) . to get 



= 1,00, = Gsjsi = V G^ - LH i , (90) 
(72 = /2^2 = Gsjc, = VG2 - , (91) 

(?3 = /3tU3 = , (92) 



insertion whereof into ( l38l) would then lead us to the same result ( l89l) . 

We would note that in the new Hamiltonian the coordinates g, h are 
cyclic, and hence the momenta G, H are integrals of motion. Also, since 

- G < L < G, (93) 
the hamiltonian is non-negative, 

7^ > (94) 

To get the canonical equations of motion, denote the generalised coor- 
dinates by q = [(?, /i, lY and the conjugate momenta by p = [G, if, L]^. 
Hamilton's equations in the absence of external torques are 

- - f (^^) 



Evaluating Hamilton's equations for (1891) . yields the canonical equations of 
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free rotational motion: 



A ^ ^^0 (98) 

ai V4 h h ) ' 

G = -^ = (100) 

L = -^ = {L^-G^)sml cosl(-~-] (102) 

Eqs. ( 1991) . (11021) are separable differential equations, and hence can be solved 
in a dosed form utilising the fact that 7i is constant (implying conserva- 
tion of energy) [2]. The {l,L) phase plane may be characterised by plotting 
isoenergetic curves of the Hamiltonian fl89l) . An example plot is depicted 
by Fig. m clearly showing the separatrix between rotational and librational 
motions. 



4 THE CANONICAL PERTURBATION 
THEORY IN APPLICATION TO 
ATTITUDE DYNAMICS AND TO 
ROTATION OF CELESTIAL BODIES 

The content of this section is based mainly on the results of papers by Efroim- 
sky [36] and Efroimsky & Escapa [37] , to which we refer the reader for details. 

4.1 A modified Andoyer set of variables 

To understand how the SA formalism may be used to model disturbing 
torques, let us start with an orbital dynamics analogy. In the theory of 
orbits, a Keplerian ellipse or hyperbola, emerging as an unperturbed two- 
body orbit, is considered as a sort of "elementary motion," so that all the 
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other available trajectories are considered to be distortions of such conies, 
distortions implemented via endowing the orbital parameters Cj with their 
own time dependence. Points of the trajectory can be donated by the "ele- 
mentary curves" either in a non-osculating manner, as in Fig. |5l3 or in the 
osculating one, as in Fig. [61 

Similarly, in attitude dynamics, a complex spin can be presented as a 
sequence of configurations constituted by some "elementary rotations." The 
easiest possibility is to use in this role the Eulerian cones, i.e., the loci of the 
spin axis, corresponding to undisturbed spin states. These are the simple 
motions exhibited by an undeformable unsupported rotator with no torques 
acting on Then, to implement a perturbed mode, we shall have to go 
from one Eulerian cone to another, just as in Figs. and E] we go from 
one Keplerian conic to another. Accordingly, a smooth "walk" over the 
instantaneous Eulerian cones may be osculating or non-osculating. 

The physical torques, the triaxiality of the rotator, and the fictitious 
torques caused by the frame noninertiality are among the possible distur- 
bances causing this "walk." The latter two types of disturbances depend not 
only on the orientation but also on the angular velocity of the body. 

In the theory of orbits, we express the Lagrangian of the reduced two- 
body problem in terms of the spherical coordinates qj = {r , (p , 6} , then 
calculate the momenta pj and the Hamiltonian 7Y(g, p) , and apply the 

^ Historically, the first attempt of using nonosculating elements dates back to Poincare 
(1897) [43j . though he never explored them from the viewpoint of a non-Lagrange con- 
straint choice. (See also Abdullah and Albouy (2001) [33], p. 430.) Parameterisation 
of nonosculation through a non-Lagrange constraint was offered in Efroimsky (2002a, b) 

[351[3Sj- 

® Here one opportunity is to use, as "elementary" motions, the non-circular Eulerian 
cones described by the actual triaxial top, when this top is unforced. Another opportunity 
is to use for this purpose the circular Eulerian cones described by a dynamically symmetri- 
cal top (and to consider its triaxiality as another perturbation) . The results of our further 
study will be independent from the choice between these two options. 
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Hamilton- Jacobi method [2H] in order to get the Delaunay constant^ 

{Qi , Q2 , Qs -Pi ) P2 , P3} = 

{ ^/Jm , \/ iia{\ — e^) , A//ia (1 — e^) cos 2 ; — Mo , — , — } , 

\i being the reduced mass. 

Very similarly, in attitude mechanics we specify a rotation state of a body 
by the three Euler angles = (p , 6 , ip and their momenta $ , , \1' . After 
that, we can perform a canonical transformation to the afore described Serret 
variables (which are, in the unperturbed case, merely constants of motion). 
A different choice of the generating function would lead one to a different set 
of arbitrary constants, one consisting of the Andoyer variables' initial values: 
{ Lo , G , H , lo , go , h} , where Lo , lo and go are the initial values of L , / 
and g . The latter selo (which we shall call "modified Andoyer set" ) consists 
only of constants of integration, and hence the corresponding Hamiltonian 
becomes nil. Therefore, these constants are the true analogues of the Delau- 
nay set with Mo (while the conventional Andoyer set is analogous to the 
Delaunay set with M used instead of Mo-). The main result obtained be- 
low for the modified Andoyer set { Lo , G , H , lo , go , h} will then be easily 
modified for the conventional Andoyer set of variables {L,G,H,l,g,h}. 

All in all, the canonical treatment of both orbital and rotational cases 
begins with 

q and p being the coordinates and their conjugated momenta, in the orbital 
case, or the Euler angles and their momenta, in the rotation case. Then one 



^ Be mindful that the more customary Delaunay variables include M instead of Mo , 
and the three latter elements have opposite signs (and play the role of coordinates, not mo- 
menta). For such variables, the Delaunay Hamiltonian is the negative of the actual Hamil- 
tonian perturbation and is, above all, amended with fj? /(2L) . However, the Hamilton- 
Jacobi procedure performed with the actual, physical Hamiltonian perturbation leads 
exactly to the Delaunay constants p03p - see section 136 in Plummer (1918). 

A similar set consisting of the initial values of Andoyer-type variables was pioneered 
by Fukushima and Ishizaki [55] . 



(103) 
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switches, by a canonical transformation 



q = /(Q, P, t) 

p = x(Q, P, t) 



(105) 



to 



« = 1)P =° . P = -aQ = « . = 

where Q and P are the Delaunay variables, in the orbital case, or the (mod- 
ified, as explained above) Andoyer variables { , C , if , /o , S'o , /i } , in 
the attitude case. 

This algorithm is based on the circumstance that an unperturbed Kepler 
orbit (and, similarly, an undisturbed Euler cone) can be fully defined by six 
parameters so that: 

1. These parameters are canonical variables { Q , P } with a zero Hamil- 
tonian: 7i*(Q, P) = ; and therefore these parameters are constants. 

2. For constant Q and P , the transformation equations (llOSp are equiva- 
lent to the equations of motion (1104p . 

4.2 The canonical treatment of perturbations 

Under perturbations, the "constants" Q , P begin to evolve, so that after 
their substitution into 

q = /(Q(t), P(t), t) 

(107) 

p = x(Q(t), P(t), t) 

(/ and X being the same functions as in f llOSp ), the resulting motion obeys 
the disturbed equations 

^ = ' P = • 
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We want our "constants" Q and P also to remain canonical and to obey 



where 

W=Q and A7i:*(Q, P t) = A7^(q(Q,P,t), p(Q,P,t), t) .(110) 

Above all, it is often desired that the perturbed "constants" Cj = Qi , Q2 , Qs 
Pi , P2 , P3 (the Delaunay constants, in the orbital case, or the modified An- 
doyer variables, in the rotation case) be osculating. This demand means that 
the perturbed velocity should be expressed by the same function of Cj{t) 
and t as the unperturbed velocity used to. In other words, the instantaneous 
"simple motions" parameterised by the "constants" should be tangent to the 
perturbed trajectory. (In the orbital case, this situation is shown on Fig. 
m) Let us check if osculation is always preserved under perturbation. The 
perturbed velocity reads 

q = g + $ (111) 

where 

„cw, , . ,,,, 

is the functional expression for the unperturbed velocity; and 

TOO -E^^^ AW (113) 

is the convective term. Since we chose the "constants" Cj to make canonical 
pairs (Q, P) obeying (11091) - (II 101) with vanishing Ti* , then insertion of 
f lTO into (flT^ will result in 

So the canonicity demand is often incompatible with osculation. Specifically, 
whenever a momentum-dependent perturbation is present, we still can use 
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the ansatz fllU7p for calculation of the coordinates and momenta, but can no 
longer use flll2p for calculating the velocities. Instead, we must use f illip . 
Application of this machinery to the case of orbital motion is depicted on 
Figini Here the constants Cj = {Qn, Pn) parameterise instantaneous el- 
lipses that, for nonzero are not tangent to the trajectory. (For details 
see Efroimsky & Goldreich [10], [H], and Efroimsky [12] •) In the case of 
rotational motion, the situation will be identical, except that, instead of in- 
stantaneous Keplerian conies, one will get instantaneous Eulerian cones (i.e., 
the loci of the rotational axis, corresponding to non-perturbed spin states). 

4.3 From the modified Andoyer variables 
to the regular ones 

Practical calculations used in the theory of planetary rotation and in space- 
craft attitude dynamics are almost always set out in terms of the regular An- 
doyer variables, not in terms of their initial values (the paper by Fukushima 
& Ishizaki [39] being a unique exception). Fortunately, all our gadgetry, 
developed above for the modified Andoyer set, stays applicable for the reg- 
ular set. To prove this, let us consider the unperturbed parameterisation 
of the Euler angles g„ = {(p , 6 , via the regular Andoyer variables 
Aj = (/, g, h- L,G,H): 

qn = fn{A,{C,t), A,{C,t)) , (115) 

each element Ai being a function of time and of the initial values Cj = 
[lo , Qo , h ] Lo , G , H ) . The explicit form of parameterisation (11151) is given 
by fl85] - [871). When a perturbation gets turned on, this parameterisation 
stays, while the time evolution of the elements A^ changes: beside the 
standard time-dependence inherent in the free-spin Andoyer variables, the 
perturbed elements acquire an extra time-dependence through the evolution 
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of their initial valuesjll] Then the time evolution of an Euler angle g„ = 
( , 6 , ip) will be given by a sum of two items: (1) the angle's unperturbed 
dependence upon time and time-dependent Andoyer variables; and (2) the 
convective term that arises from a perturbation-caused alteration of the 
Andoyer variables' dependence upon the time: 

g„ = g„ + . (116) 

The unperturbed part is 

while the convective term is given by 



where the set Cj is split into canonical coordinates and momenta like this: 
Qj = y do, h) and Pj = {Lo , G , H) . In the case of free spin they 
obey the Hamilton equations with a vanishing Hamiltonian and, therefore, 
are all constants. In the case of disturbed spin, their evolution is governed by 
(11091 - [TTO|) . substitution whereof in (11181) once again takes us to (11141) . This 
means that the nonosculation-caused convective corrections to the velocities 
stay the same, no matter whether we parameterise the Euler angles through 



This is fuUy analogous to the transition from the unperturbed mean longitude, 
M (t) — Mo + n (t — to) , with Mq , n , to — const , 
to the perturbed one, 

M{t) = Mo{t) + n{t') dt' , with to = const , 

in orbital dynamics. 
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the modified Andoyer variables (variable constants) or through the regular 
Andoyer variables. This invariance will become obvious if we consider the 
analogy with orbital mechanics: in Fig. [5], the correction $ is independent 
of how we choose to parameterise the non-osculating instantaneous ellipse 
- through the set of Delaunay constants containing Mo , or through the 
customary set of Delaunay elements containing M . (The latter set is, 
historically, chosen to correspond to a Hamiltonian variation taken with an 
opposite sign - see footnote 8 above. This subtlety, however, is unimportant 
to our point.) 



4.4 The Andoyer variables introduced 
in a precessing frame of reference 

4.4.1 Physical motivation 

Let us consider when the perturbing torque depends not only on the 

instantaneous orientation but also on the instantaneous angular velocity of 
the rotator. In particular, we shall be interested in the fictitious torque 
emerging when the description is carried out in a precessing coordinate sys- 
tem. This situation is often encountered in the theory of planetary rotation, 
where one has to describe a planet's spin not in inertial axes but relative to 
a frame associated with the planet's circumsolar orbital plane (the planet's 
ecliptic). The latter frame is noninertial, because the ecliptic is always pre- 
cessing due to the perturbations exerted by the other planets. The reason 
why astronomers need to describe the planet's rotation not in an inertial 
frame but in the ecliptic one is that this description provides the history of 
the planet's obliquity, i.e., of the equator's inclination on the ecliptic. As the 
obliquity determines the latitudinal distribution of insolation, the long-term 
history of the obliquity is a key to understanding climate evolution. Interest- 
ingly, the climate is much more sensitive to the obliquity of the planet than 
to the eccentricity of its orbit. (Murray et al Ward |18], [19]) 

The canonical theory of rotation of a rigid body in a precessing coordi- 
nate frame was pioneered by Giacaglia and Jefferys [31] . It was based on the 
Andoyer variables and was used by the authors to describe rotation of a space 
station. This theory was greatly furthered by Kinoshita [8] who applied it 
to rotation of the rigid Earth. Later it was extended by Getino, Ferrandiz 
and Escapa [9], [10], [11] to the case of nonrigid Earth. Simplified versions 
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of the Kinoshita theory were employed by Laskar and Robutel [H] and by 
Touma and Wisdom [TB] in their studies of the long-term evolution of 
the obliquity of Mars. While a detailed explanation of this line of research 
will require a separate review paper, here we shall very briefly describe the 
use of the Andoyer variables in the Kinoshita theory, and shall dwell, follow- 
ing Efroimsky & Escapa on the consequences of these variables being 
nonosculating. 

4.4.2 Formalism 

Consider an unsupported rigid body whose spin should be described with 
respect to a coordinate system, which itself is precessing relative to some 
inertial axes. The said system is assumed to /X so the 

kinetic energy of rotation, in the inertial axes, is given by 



where uj is the body-frame-related angular velocity of the body relative to 
the precessing coordinate system, while /x is the angular velocity relative 
to the inertial frame. In ( 1119^ . both uj and /j, are resolved into their 
components along the principal axes. The role of canonical coordinates will 
be played the Euler angles = 4> , G , that define the orientation of the 
principal body basis relative to the precessing coordinate basis. To compute 
their conjugate momenta p„ = $ , 6 , \1/ , let us assume that noninertiality 
of the precessing coordinate system is the only angular-velo city- dependent 
perturbation. Then the momenta are simply the derivatives of the kinetic 
energy. With aid of ([S]-[7]), they can be written as 

$ = ^ = Ji {uJi + /ii) sin 6' siiaip + I2 (tU2 + /^2) sin 6' cosip + I3 (tus + /is) cos 6' , (120) 




(119) 



d(f) 



e 



de 



h (t^3 + yWs) 



(121) 



dip 



h {oJi + /ii) cosip - I2 {uj2 + /i2) sin-?/' 



(122) 
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These formulae enable one to express the angular-velocity components oOi 
and the derivatives = {(p , 6 , ip) via the momenta = ( $ , , \1/ ) . 
Insertion of ([I20]-[I^ into 

^ = 5^ qnPn - £ = 0$ + + - T + 9, ij) (123) 



results, after some algebra, in 

H = T + AH 

where 



(124) 



AH 



/ii 



/i2 



sin ^ 

cosip 
sin ^ 



($ - ^ cos6') + e cosV' 
($ - \1/ cos 61) - 6 sin?/' 



(125) 



and the potential V is presumed to depend only upon the angular coordi- 
nates, not upon the momenta. As can be easily seen from fl25] - [27|) . formula 
(11251) is simply another form of the relation ATi. = — fi - g which, according 
to (jnni - E^]) , can also be expressed via the Andoyer variables solely: 



AH 



/i 



sin/ — fi2 y/G'^ — cos / — fi^ L . (126) 



Now let us employ the machinery described in the preceding subsection. 
The Euler angles connecting the body axes with the precessing frame will 
now be expressed via the Andoyer variables by means of flll5p . (The explicit 
form of the functional dependence (11151) is given by (IHSl- IHTl) . but this exact 
form is irrelevant to us.) The fact that the Andoyer variables are introduced 
in a noninertial frame is accounted for by the emergence of the /x-terms in 
the expression (11251) for the disturbance ATi. . Insertion of (11251) into (IllSp 
entails: 



Qn = g„ + 



dAH 



the convective terms being given by 

dAH dAH /ii simp + ^2 cosip 



dPi 



9$ 



sin 6' 



(127) 



128) 
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dAH dAH , . , 

"qP^ = ~qq' = ~ /^i co^V^ + simjj , (129) 

dAH dAH , . , ,N /.onN 

~apr " " [f^i smtp + H2 cosf) cote - fi3 . (130) 

It should be stressed once again that the indices n = 1 , 2 , 3 in (11271) 
number the Euler angles, so that (jn stand for (p , 9 , , and pn signify 
$ , 9 , \E' . At the same time, the subscripts i = 1,2,3 accompanying the 
components of /J, in (11281) correspond to the principal body axes. 

4.4.3 The physical interpretation of the Andoyer variables defined 
in a precessing frame. 

The physical content of the Andoyer construction built in inertial axes is 
transparent: by definition, the element G is the magnitude of the angular- 
momentum vector, L is the projection of the angular-momentum vector on 
the principal axis ba of the body, while H is the projection of the angular- 
momentum vector on the §3 axis of the inertial coordinate system. The 
variable h conjugate to H is the angle from the inertial reference longitude 
to the ascending node of the invariable plane (the one perpendicular to the 
angular momentum). The variable g conjugate to G is the angle from the 
ascending node of the invariable plane on the reference plane to the ascending 
node of the equator on the invariable plane. Finally, the variable conjugate 
to L is the angle / from the ascending node of the equator on the invariable 
plane to the the bi body axis. Two auxiliary quantities defined through 

H L 

too, have evident physical meaning: / is the angle between the angular- 
momentum vector and the §3 space axis, while J is the angle between the 
angular-momentum vector and the b3 principal axis of the body, as can be 
seen on Fig. [2J 

Will all the Andoyer variables and the auxiliary angles I and J retain 
the same physical meaning if we re-introduce the Andoyer construction in 
a noninertial frame? The answer is affirmative, because a transition to a 
noninertial frame is no different from any other perturbation: precession of 
the fiducial frame (§1 , §2 , §3) is equivalent to emergence of an extra per- 
turbing torque, one generated by the inertial forces. In the original Andoyer 
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construction assembled in an inertial space, the invariable plane was orthog- 
onal to the instantaneous direction of the angular-momentum vector: If the 
perturbing torques were to instantaneously vanish, the angular-momentum 
vector (and the invariable plane orthogonal thereto) would freeze in their 
positions relative to the fiducial axes (si , §2 , §3) (which were inertial and 
therefore indifferent to vanishing of the perturbation). Now, that the An- 
doyer construction is built in a precessing frame, the fiducial plane is no 
longer inertial. Nevertheless if the inertial torques were to instantaneously 
vanish, then the invariable plane would still freeze relative to the fiducial 
plane (because the fiducial plane would cease its precession). Therefore, all 
the variables will retain their initial meaning. In particular, the variables 
/ and J defined as above will be the angles that the angular-momentum 
makes with the precessing §3 space axis and with the b3 principal axis of 
the body, correspondingly. 

4.4.4 Calculation of the angular velocities via the Andoyer 
variables introduced in a precessing frame of reference 

Let us now get back to formulae (0 - E]) for the principal-body-axes-related 
components of the angular velocity. These formulae give the angular velocity 
as a function of the rates of Euler angle's evolution, so one can symbolically 
denote the functional dependence ([5] -[7]) as uj = t<;(g) . This dependence is 
linear, so 

u;{qiA)) = a;(g(A)) + cjidAH/dp) , (131) 

A being the set of Andoyer variables. Direct substitution of (11281 - [T30|) into 
((5] - [7]) will then show that the second term on the right-hand side in (11311) 
is exactly — /j. : 

u:{q{A)) = u;{g{A)) - /x • (132) 

Since the "total" angular velocity uj{q) is that of the body frame relative to 
the precessing frame, and since /i is that of the precessing frame relative to 
some inertial frame, then ijj{g{A) ) will always return the angular velocity 
of the body relative to the inertial frame of reference, despite the fact that the 
Andoyer variables A were introduced in a precessing frame. This nontrivial 
fact has immediate ramifications for the theory of planetary rotation. These 
will be considered in the subsequent subsections. 
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In short, the above may be summarised as: 
= a;(g(A)) + u:{d/S.n/dp) , 

uj{q{A)) = ^'"^ , 
u){dAH/dp) = At , 



a; 



(re!) 



(inert) 



where the entities are defined as follows: 



(rel) 



uj{g{A)) = UJ 



{inert 



UJ = the relative angular velocity, 

i.e., the body's angular velocity relative to a precessing orbital frame 



the precession rate of that frame with respect to an inertial one; 



UJ 



(inert) 



the inertial angular velocity, 

i.e., the body's angular velocity with respect to the inertial frame. 



To better understand the origin of fll33p . let us get back to the basic Andoyer 
formalism introduced in the previous subsection. As the first step, we intro- 
duce, in an unperturbed setting (i.e., for an unsupported rigid rotator con- 
sidered in an inertial frame), the parameterisation of the Euler angles g„ = 
( , 9 , ip) through the Andoyer variables Aj = {I , g , h; L , G , H) : 

qn = fn{A,{C,t), A,{C,t)) , (134) 

each variable Ai being dependent upon its initial value Ci and the time: 



Ai — d 



rij dt 



(135) 



to 



where the mean motion^ rii may bear dependence upon the other Andoyer 
variables Aj . As already mentioned above, this is analogous to the evolution 
of the mean longitude M = Mq + n{t — to) in the undisturbed two-body 



Three of the six Andoyer variables have vanishing Ui in the unperturbed free-spin 
case, but at this point it is irrelevant. 
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problem (where the mean motion n is a function of another orbital element, 
the semimajor axis a ). 

The initial values Ci are integration constants, so in the unperturbed 
case one can calculate the velocities Qn = (j) , 9 , ipj simply as the partial 
derivatives 

As the second step, we employ this scheme in a perturbed setting. In 
particular, we introduce a perturbation caused by our transition to a coor- 
dinate system precessing at a rate . Our Euler angles now describe the 
body orientation relative to this precessing frame. By preserving the pa- 
rameterisation (I134p . we now introduce the Andoyer variables Ai in this 
precessing frame. Naturally, the time evolution of Ai changes because the 
frame-precession-caused Hamiltonian disturbance ATi now shows itself in 
the equations of motion. This disturbance depends not only on the body's 
orientation but also on its angular velocity, and therefore our Andoyer vari- 
ables cannot be osculating (see equation flll4p and the paragraph thereafter). 
This means that the unperturbed velocities, i.e., the partial derivatives (I136p 
no longer return the body's angular velocity relative to the precessing frame 
(i.e., relative to the frame wherein the Andoyer variables were introduced). 
This angular velocity is rather given by the sum flll6p . However, as explained 
above, the unperturbed expression g{A) , too, has a certain physical mean- 
ing: when plugged into <-i;(g) , it always returns the angular velocity in the 
inertial frame. It does so even despite that now the Andoyer parameterisa- 
tion is introduced in a precessing coordinate frame 1^ 

This parallels a situation in orbital dynamics, where the role of canonical variables 
is played by the Delaunay constants - see expression (|103p and a footnote accompanying 
it. In the unperturbed setting (the two-body problem in inertial axes), the Cartesian 
coordinates r = ( xi , a;2 , xa) and velocities (ii , ±2 , is) are expressed via the time 
and the Delaunay constants by means of the following functional dependencies: 

^ = f {C , t) and V = g ( C , t ) , where g = df /dt . 

If we want to describe a satellite orbiting a precessing oblate planet, we may fix our 
reference frame on the precessing equator of date. Then the two-body problem will get 
amended with two disturbances. One, ATiotiate , caused by the presence of the equa- 
torial bulge of the planet, will depend only upon the satellite's position. Another one, 
ATipj-ecess 1 will stcm from the noninertial nature of our frame and, thus, will give birth 
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4.4.5 Example 1. The theory of Earth rotation 

As an example, let us consider the rigid-Earth-rotation theory by Kinoshita 
[8]. Kinoshita began with the standard Andoyer formalism in inertial axes. 
He explicitly wrote down the expressions flllSp for the Euler angles g„ = 
fn{A) of the figure axis of the Earth, differentiated them to get the ex- 
pressions for the velocities g„ as functions of the Andoyer variables: (jn = 
gn{A) , and then used those expressions to write down the angles Ir{g{A)) 
and hr{g{A)) that define the orientation of the rotation axis0 Kinoshita 
pointed out that one's knowledge of the solar and lunar torques, exerted 
on the Earth due to its nonsphericity, would enable one to write down the 
appropriate Hamiltonian perturbations and to calculate, by the Hori-Deprit 
method ([3] - [7]), the time evolution of the Andoyer variables. Substitution 
thereof into (11341) would then give the time evolution of the Euler angles 
that define the figure axis of the planet. Similarly, substitution of the cal- 
culated time dependencies of the Andoyer variables A into the expressions 
for Ir{g{A)) and hr{g{A)) would yield the time evolution of the planet's 
rotation axis. 

The situation, however, was more complicated because Kinoshita's goal 
was to calculate the dynamics relative to the precessing ecliptic plane. To 
achieve the goal, Kinoshita amended the afore described method by adding, 
to the lunisolar perturbations, the momentum-dependent frame-precession- 
caused term (see our formula (11251) above). Stated differently, he introduced 
the Andoyer variables in a precessing frame of reference. This made his An- 
doyer variables A nonosculating. Kinoshita missed this circumstance and 
went on to calculate the time dependence of the so introduced Andoyer vari- 
ables. He then plugged those into the expressions = fn{A) for the Euler 

to velocity-dependent inertial forces. Under these perturbations, the Delaunay constants 
will become canonical variables evolving in time. As explained in subsection 4.2, the 
velocity-dependence of one of the perturbations involved will make the Delaunay variables 
nonosculating. On the one hand, the expression r = f (C{t) , t) will return the cor- 
rect Cartesian coordinates of the satellite in the precessing equatorial frame, i.e., in the 
frame wherein the Delaunay variables were introduced. On the other hand, the expression 
g ( C , t) will no longer return the correct velocities in that frame. Indeed, according to 

- I114p , the Cartesian components of the velocity in the precessing equatorial frame 
will be given by g(C, t) + dAHprecess/dp . However, it turns out that g(C, t) 
renders the velocity with respect to the inertial frame of reference. (Efroimsky 42], [50] ) 

For the Euler angles, Kinoshita chose notations, which are often used by astronomers: 
h , I , <j) and which are different from the convention = (j) ^ 9 , ip used in physics. 
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angles of the figure and into the expressions Ir{g{A)) and hr{g{A)) for the 
orientation angles of the Earth rotation axis. The expressions g„ = fn{A) 
still gave him the correct Euler angles of the Earth figure (now, relative to 
the precessing ecliptic plane). The expressions Ir(g{A)) and hr{g{A)) did 
NOT give him the correct orientation of the angular-velocity vector relative 
to the precessing frame, because the Andoyer variables introduced by Ki- 
noshita in the precessing frame were nonosculating. This peculiar feature of 
his theory has long been ignored, because no direct methods of measurement 
of the angular velocity of the Earth had been developed until 2004. While the 
thitherto available observations referred only to the orientation of the Earth 
figure (Kinoshita [51]), a technique based on ring laser gyroscopes (Schreiber 
et al. [52]) and a VLBI technique (Petrov 2007 [53]) later made it possible 
to directly measure the instantaneous angular velocity of the Earth relative 
to an inertial frame. 

As we demonstrated above, when the Andoyer variables A are intro- 
duced in a precessing frame, the expressions = fn{,A) return the Euler 
angles of the body relative to this precessing frame, while uj{g^{A) ) re- 
turns the body-frame-related angular velocity relative to the inertial frame. 
Accordingly, the expressions Ir{g{A)) and hr{g{A)) return the direction 
angles (as seen by an observer located on the body) of the instantaneous 
angular velocity relative to the inertial frame, i.e., the velocity observed in 
[52]. This way, what might have been a problem of the canonical Kinoshita 
theory became its advantage. It is exactly due to the nonosculation of the 
Andoyer elements, introduced in a precessing ecliptic frame, that this theory 
always returns the angular velocity of the Earth relative to inertial axes. We 
see that sometimes loss of osculation may be an advantage of the theory. 

4.4.6 Example 2. The theory of Mars rotation. 

Two groups ([H], [151 [IS]) independently investigated the long-term evolu- 
tion of Mars' rotation, using a simplified and averaged version of the Ki- 
noshita Hamiltonian. The goal was to obtain a long-term history of the 
Martian spin axis' obliquity, i.e., of the angle between the Martian spin axis 
and a normal to the Martian ecliptic. Just as in the afore described case of 
the Earth, the Martian spin axis is evolving due to the solar torque acting 
on the oblate Mars, while the Martian ecliptic plane is in precession due to 
the perturbations exerted upon Mars by the other planets. Since for realistic 
rotators the Andoyer angle J is typically very small (i.e., since the angular- 
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velocity and angular-momentum vectors are almost parallel), one may, in as- 
tronomical applications, approximate the obliquity with the angle made by 
the planet's angular-momentum vector and the ecliptic. When the Andoyer 
construction is built in a precessing frame (the fiducial basis (§1 , §2 , §3) 
being fixed on the planet's ecliptic), this assertion means that the obliquity 
is approximated with the Andoyer angle / . When the Andoyer variables 
are introduced in the traditional way (the basis (§1 , §2 , §3) being inertial), 
the above assertion means that one has to find the orientation of the angular 
momentum relative to the the inertial axes (i.e., to calculate the Euler angles 
h , I ) and to find the orientation of the ecliptic relative to the same inertial 
axes. Then the orientation of the angular momentum with respect to the 
ecliptic will be found, and it will be an approximation for the obliquity. 

The former approach was implemented by Laskar and Robutel [l^, the 
latter by Touma and Wisdom [151 EH] . Despite that these two teams intro- 
duced the Andoyer variables in different frames of reference, the outcomes 
of their calculations were very close, minor differences being attributed to 
other reasons!^ This coincidence stems from the afore explained fact that 
the Andoyer variables and the Andoyer angles / , J retain their physical 
meaning when introduced in a precessing frame. (See subsection 4.4.3.) 

5 The Sadov Variables 

The Andoyer variables have the merit that they reduce the Hamiltonian of an 
unsupported and torque-free rigid body to one and a half degrees of freedom, 
the resulting expression for the Hamiltonian being very simple. However, 
except for the case of axial symmetry, these variables are not action-angle. 

Sadov [22] and Kinoshita [27j, in an independent way, obtained sets of 
action-angle variables for the rotational motion of a triaxial rigid body. Both 
sets of variables are very similar. Essentially, they are obtained by solving the 
Hamilton- Jacobi equation stemming from the Hamiltonian in the Andoyer 
variables. Sadov's transformation is formulated in terms of the Legendre 
elliptic functions of the first and third kind, while in the Kinoshita transfor- 
mation the Heuman Lambda function emerges. 

Here we present the approach of Sadov who began by introducing an 
intermediate set of canonical variables (/? , a) that nullify the Hamiltonian. 

While Touma and Wisdom [TSKTB] employed unaveraged equations of motion, Laskar 
and Robutel [13] used orbit-averaged equations. 
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The generating function S of the transformation 

{£, g, h, L, G, H) — > {(3i, (32, Ps, oi, 02, "3), 
can be found by solving the corresponding Jacobi equation 



1 f sinH cos^£ 



/dS_ 



as 



1 fdsy dS ^ 



Since the system is autonomous and variables g and h are cyclic, the gener- 
ating function may be expressed as 



S — —ait + a2g + a^h + U{£;ai,a2,as) , 



wherefrom 



dU 



de 



1 

^3 



1 

T2 



1 

h 



2a\ - a\ — 



-an 



and further simplification yields a simple equation for the function lA : 

'dlj ~ ^ [c + dsinH 



with 



a = h{al - 2ail2) , c = h{h - h) > 



b={l2- h)al > , hih - /i) > . 

Hence, the transformation is 



a + b sin^ i 
c + dsm^e ' 



dU 
dai ' 



OS 

G = — = a2 , 
og 



p2 = g + 



dU 
da2 



(137) 



whence we see that 
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ai. 



(138) 



At this point, Sadov introduceco a parameter k and a so-called state 
function A as 



«:^ = ^4^^>0, X^ = n^'-1^4^1^>0 . (139) 



Via these quantities, the angular momentum L in (11371) may be expressed as 

'2 /c2 + A2 1 + /t2 sin^ 



/.^ l-A^ + (/.^ + A^)sin^£ 



and the Hamiltonian (I138p becomes: 

2/1/3 «:' + A2 ■ 

With this transformation accomplished, Sadov proceeded to a second one, 

(£, 5(, /i, L, G, //) — ^ {ipe,, (fg, (fh, Ig, h) , (141) 
the new actions being 

L = — (D Ldi, la = — (p G dq = G = ao, h = — (p H dh = H = a^i 
27r / ' ^ 271 J ^ ' 2n J 

While the variables (/3 , a) corresponded to a vanishing Hamiltonian, the 
action-angle ones correspond to the initial Hamiltonian of Andoyer (though 
now this Hamiltonian has, of course, to be expressed through these new 
variables) . 

By means of a convenient auxiliary variable z introduced through 



cos 2; „ vT+K^sinz , , , 

sin£ = ^^ , cos£ = -^ , (142) 

V 1 + sin z V 1 + ^2 sin z 

it is possible to derive from (I140p that 



Ie=^hdi= ^^^f^ [{k' + A^)n(.^ A) - A^ K{X)] . (143) 

27r J TTKy/K^ + 

where K and H are the complete elliptical integrals of the first and the third 
kind, respectively. 



16 



These quantities already appear, in a different context, in f25', p. 394 
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Symbolically, the above expression may be written as = Ie{Ig , A) = 
Ig /(A) . According to the Implicit Function Theorem, there exists a func- 
tion (f){l£/Ig) inverse, locally, to /(A). Let us denote this function as 
A = (f){Ie/Ig) . Sadov ^29j proved that it is defined for all values {h/ Ig) > , 
and that it is analytic for [h/ Ig) > or {h/ Ig) 7^ (S/vr) arctauK at A = 1 . 

Since the Hamiltonian is known to be 

the angular variables could be found by integrating the Hamilton equations 

du on dn 

if the explicit relation between the action momenta and A were known. Un- 
fortunately we do not have such a relation at hand. To circumvent this 
difficulty, the generating function of the transformation must be derived. 

Let W be the generating function of the transformation (11411) from the 
Andoyer variables to the action and angle variables. It may be chosen as 

W = Igg + hh + V{t,hJg) , 

and the equations of this transformation will read: 



dV 

dig 

h . 

The generating function V is obtained directly from the quadrature 

V{t,h,Ig) = f Li£;k,Ig)di , 
Jio 

where we have to replace L with expression fll40p . 

To calculate ipe and ipg , we need to find the partial derivatives dX/dli 
and dX/dIg . Since A = (j){Ie/Ig) and le = Igf{X) , then 

dh _ df _ 1 rdL_ /,2/tA(l + ^^)^/^ 
dX 'dX 2nJdX 7r(/t2 + A2)3/2 





dV 






dV 


L = 


'de 




V^e = 


We 


G = 








dW 


dg 




Vg = 


dig 


H = 


dW 






dW 


dh 


= 4 , 


Vg = 


dh 
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After introducing notation u = h/ Ig , we can write: 

d(P _ 1 _ 1 _ 1 rrjK^ + Xy/^ 

~ du/d(P ~ ^dle/d(f) ~ ^dh/dX ~ ~ 2kX{1 + k^Y/^K{\) 



(145) 



whence 
dX_ _ c 



du 



+ X^ 



du dig n 2kA(1 + K^y/^K{X) IgK^XK{X) 



[{k^ + X^)U{k\X)-X^K 



Now, taking into account fll4UI) and employing the chain rule, we arrive at 



dV dV dX dX 



dip dX dl 



dh 



dL 



^0 



dX 



¥^g= 9 + 



_d_ 
dl 



Ldi = g + 



9 Jio 



Ldi + 



'9 Je.0 



dX_ 
dl 



9 J to 



dL 

'dX 



di . 



Straightforward differentiation of (11401) yields: 
dL _ IgKX{l + K^) 1 



dX 



(^2 + ^2)3/2 Vl + ^2 sin^ iJl - A2 + (k^ + X^) siu" 



With the help of (I142p . we can now compute the quadratures 

' dL IgKXil + K^fl^ 



dX 



(^2 + ^2)3/2 



F(z,A) 



and 



Ldl 



l + /€^ 



insertion whereof into Eqs. (I146p results in 
7rF(z,A) 



/?2 + A')n(2, h',X)-X^K{X) 



2 K{\) ' 

g+-^{K^ + X^){l + K^) 



K 



U{Z,K^,X) - 



U{k^,X)F{z,X) 



K{\) 



Inversion of the first of these expressions entails 

z = am{2{piK{X)/n, A) . 



(146) 



(147) 
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(For the definition of tlie elliptic function am see Appendix A. 2 below. ) 
Thus, we have the Andoyer angle g expressed via the action-angle variables: 

while the angle i is obtained directly from expression fll42p . 

6 CONCLUSIONS 

In this paper, we have reviewed the Serret- Andoyer (SA) formalism for mod- 
eling and control of rigid-body dynamics from the dynamical systems per- 
spective. We have dwelt upon the important topic of modeling of general, 
possibly angular velocity-dependant disturbing torques, and upon the inter- 
connection between the Andoyer and the Sadov sets of variables. We have 
also contributed some new insights. 

The first insight is that the Andoyer variables turn out to be non-osculating 
in the general case of angular-velocity-dependent perturbation. The second 
insight is that even when these variables are introduced in a precessing ref- 
erence frame, they preserve their interconnection with the components of 
the angular momentum - a circumstance that makes the Andoyer variables 
especially valuable in astronomical calculations. 

In summary, this treatise constitutes a first step towards understanding 
the consequences of using the SA formalism as a single, generic language 
for modelling rigid body dynamics in diverse - and seemingly unrelated - 
areas such as celestial mechanics, satellite attitude control, and geometric 
mechanics. 

Appendix. 



A.l. Spherical-trigonometry formula (1571) 

The standard formula of spherical trigonometry |34j . 

cosg = cos(0 — h) cos{tp — + sin(0 — h) sm{ip — I) cos{n — 6) , (148) 



n( am(2^^ir(A)/7r, A),^^, A 



42 



immediately entails: 

singdg = [sin(0 — /i) cos^ip — I) — sin(0 — /i) sin(-?/' — /) cos 6 ] d{(j) — h) 
+ [cos(0 — /i) sin{ip ~ I) — sin(0 — /i) cos{tp — I) cos 9 ] d{4' — I) 

+ sin(0 — /) sin{ilj — I) sin6 dO . 

The other standard formulae of the spherical trigonometry enable one to 
transform (11491) into 

sin g dg = 

smg cosId{(j) — h) + sing cosJ d^ip — l) + sui{(j) — h) sin^ip — l) sin 9 dO 

sing cosId{(t> — h) + sing cosJd{ip — l) + sin{i/j — I) sin J sing dO , 
wherefrom equality (1571) follows. 

A. 2. The dimensionless time u and its interconnection 
with Tkin , G , G cos / , 9 , and i/j 

First let us recall some basics. As agreed in the text, we choose the body 
axes to coincide with the principal axes of inertia, ti , la2 , bs . This makes 
the inertia tensor look like: 



Decomposing the body angular velocity oj over the principal basis. 




Ji 
/2 
/3 



(151) 



a; = tui bi + co'2 b2 + cus bs 



(152) 



we can write twice the kinetic energy as 



2Tr 



kin 



Iiiul + I2UJ2 + hojl 



(153) 
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and the body-frame-related angular-momentum vector as 

g = I-u} . (154) 

The square of this vector will be: 

g2 = If oof + Ifuji + Ilool , (155) 

while its direction cosines with respect to an invariable plane-based coordi- 
nate system will read0 

, gl h ^1 o' §2 h ^2 ' gs ^3 

« = 77 = —PT- , ^ = r = ' 7 = 77 = , (156) 

G denoting the magnitude of the angular momentum vector: G = |g| . At 
this point it is convenient to define an auxiliary quantity via 

h h I3 

Substitution of 01561) into 01571) shows that this auxiliary quantity obeys 

p2 = h oof + hoof + I3 uol (158) 

or, equivalently, 

G^P^ = 2Tkin . (159) 

We see that P^ is an integral of motion, a circumstance that will later help 
us with reduction of the problem. 
Due to the evident identity 

{a'f + iP'Y + (7? = 1 (160) 



The direction cosines can be expressed in terms of the angle I from the ascending 
node of the equator on the invariable plane to the the bi body axis, and the angle J 
between the angular- momentum vector and the bs principal axis of the body (Fig.[2|): 

a = sin J sin I , P = sin J cos / 7 = cos J , 

so that 

IiLUi — G sin J sin I , I2UJ2 = G sin J cos I , I^oj^ = G cos J . 



44 



our depends upon only two directional cosines. Elimination of 7' from 
(11571) . by means of fll60p . trivially yields: 



Note that, if we now introduce an and a 6^ such that 

= n > (162) 



(163) 



we shall be able to write down the definition of in a form that will, 
formally, be identical to definition of an ellipse: 

+^ = 1 . (164) 



p2 _ 


1 


1 


1 


h 


h 




1 

T, 


1 


1 


h 


h 



a2 62 

This means that, along with P^ , it is convenient to introduce another 
auxiliary variable, ^ , one that obeys 

= cos^e , ^ = sin^e . (165) 



a 



2 



Insertion of the above into (I160p entails 

(7? = (1 - a^)(l - ^^si^'^) • (166) 

Then, after defining the quantity and the function A as 

6^ - 

. (167) 

and 



= - n"^ sin^^ , (168) 
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we shall be able to cast the direction cosines of the angular-momentum vector 
in the form of 



a' = a cos^ , (169) 
(3' = hsmi , (170) 

7' = Vl . (171) 

Looking back at (I156p . we see that initially we started out with three direction 
cosines, only two of which were independent due to the equality ( I160p . The 
latter meant that all these cosines might be expressed via two independent 
variables. The quantities and were cast for the part. (Mind that a 
and h depend on through (11621) - (11631) .) It should also be mentioned 
that, due to (I160p . P^ is a constant of motion, and therefore formulae (11691 
- I17ip effectively reduce the problem to one variable, ^ , which thereby plays 
the role of re-scaled time, in terms whereof the problem is fully solved. Below 
we shall explicitly write down the dependence of the "time" ^ upon the real 
time t or, equivalently, upon the dimensionless time u = n(t — to) 
emerging in (161) - [62]) . 

Undisturbed spin of an unsupported rigid body obeys the Euler equations 
for the principal-axes-related components of the body angular velocity: 

I,^ = uj, uj, {h - h) , (172) 



I,^ = 00, (J3 - h) , (173) 
dt 



I,^ = U2 (/i - I2) . (174) 
dt 

Insertion of expressions (11561) into (11721) . (11731) . and (11741) provides an equiv- 
alent description written in terms of the angular-momentum vector's direc- 
tional cosines relative to the principal axes: 

G & 1 - ^\ , (175) 



dt \h h 
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'2 

Substitution of ([163 -[HI]) into (ITTHD will then entail 



A subsequent insertion of (11621 - [T63l) and (11681) into (I178P will then entail: 



a/1 - ^2 sin^,^ V V-^i / ^3 

Now define a "mean motion" n as 



Since, according to (11591) . is integral of motion, then so is n , and 

therefore 



/ /, , • 2 =7 = / ridt' = n{t - to) 
Jo \/\ — sm 4 

This means that, if we define a dimensionless time as 



;i8i) 



u = I n dt' = n {t - to) (182) 

J to 



and a function F as 

^0 a/I - ^2 sm''^' 
then the interrelation between the parameter ^ and the time will look like: 

u = F{^, k) . (184) 
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This is a Jacobi elliptic equation whose solutions are be written in terms of 



the following elliptic functions: 

^ = am(M, k) = F-\u, k) , (185) 

sin^ = sn('U , k) = sin(am(M, k) ) , (186) 

cos^ = cn(-u , k) = cos(am(M, k) ) , (187) 

= V^l - ^2 sn2(M, k) = dn{u, k) . (188) 



According to (11841) . u is a function of ^ and k. The integral of motion k 
is, through (11671) and (11591) . a function of Tkm and G , while the parameter 
C, is, through the medium of (11621 - [T63|) and (11691 - [TTOl) . a function of Tkin , 
G , a' = sin J sin I , and (3 = sin J cos / . All in all, m is a function of 
Tkin , G , J , and / . 

It is also possible to express u through another set of geometric vari- 
ables. Recall that ( , 6 , ip) are the Euler angles defining orientation of 
the body axes ( bi , b2 , los ) relative to the fiducial frame ( §1 , §2 , §3 ) , 
while ( 00 , J , Euler angles defining orientation of the body 

axes with respect to the invariable plane, as on Fig. 2. Out of that picture, 
it is convenient to single out a spherical triangle whose sides are given by 
(0 — h , ip — I , g) , and whose internal angles are ( / , vr — 9 , J) , as 
shown in Fig. 3. An analogue to the law of cosines, for spherical triangles, 
looks like: 

cos/ = cos 9 cos J + sin 9 sin J cos('0 — /) , (189) 

which is the same as (11481) . up to a cyclic transposition. The second term on 
the right-hand side of (11891) can be expanded as 

sin J cos^ip — I) = sin J {cosip cos/ + sin'0 sin/) 

= ( sin J sin / ) sin ip + ( sin J cos I ) cos ip 
= a' sin -0 + /? cos ip 

= a sin'0 cos( am(M)) -|- b costp sm{am{u)) , 

(190) 
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where the last hne was obtained with aid of 01691 - I170p and fll85p . Thus 
(11891) acquires the shape of 

cos/ = cos^ Vl — Aam(M) + a sin-?/' cos(am(M)) 

+ h cosip sin(am(M)) . (191) 

We see that u depends upon ip , 9 , cos I , and (through a and b ) 
upon T^in and G . This is equivalent to saying that m is a function of 
Tkin , G , G COS / , 9, ^ . 



A. 3. Taking variations of S 

Let us start with the expression (158|) for the generating function: 

S = - tTkin + G h cos/ + G g + G j cos. J dl + C . (192) 

We see that it depends upon eight variables, some of which are dependent 
upon others. In brief, 

S = S{t , Tkin , I , L = G cos J , g , G , h , H = G cos / ) . (193) 
Variation thereof, taken at a fixed time t , will read: 

= ( Sn^n +(^)sh+( ..TT^^) 6{G COS/) 



dTkrnJ \dhj \d{G cos/) 



— t STkin + G cos / 6h + h 6{G cos/) 



+ [g + I cosJdl]6G + G6g + G I 5(cos J) dl . (194) 
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Our goal is to simplify fll94p . having in mind that the variables emerging 
there are not all mutually independent. To that end, let us employ expression 
(11501) . Its differentiation would yield the expression 

d g — d {4>~h) cos I = — cos J d{l — 11))+ sin J sin (/ — ip) dO , (195) 

where the Euler angles 0, 6, ip determine the orientation of the body relative 
to some inertial reference frame. If, however, we perform variation of (I150p . 
for a fixed orientation of the body relative to the inertial frame, then we shall 
get simply 6g = — cos J 61 — cos / 6h . Multiplying this by G, we arrive at 

G 6g = - G cos J 61 - G cosi 6h , (196) 

substitution whereof into fll94p entails 

6S = -t6Tkin + h6{GcosI) + (^g + jcosJdl^ 6G 

+ G j 6{cosJ) dl - G cos J 61 . (197) 

Next, we shall recall the formula for the Binet Ellipsoid, expressed through 
the directional cosines defined in fll59[ According to fll57p and (11590 . we 

Insertion of the (evident from Fig. [2]) relations 

a' = sin J sin / , = sin J cos I 7 = cos J , 

into the above formula will yield 

2 J ( sin^ / cos^ / \ cos^ J 2 Tfcj„ 



or, equivalently, 



2x I sin^l cos^l l\ 2 Tkin 1 



:i-WY){—^^~-)^^-- . (200) 
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Differentiation of the above will lead us to 



„ dTkin . Thin dG „ ' ,/ 'n / sin^ Z cos^ Z 1 ' 



+ 2 a0 \-^-]-]dl . (201) 



h I 



If we now turn the differentials into variations, and multiply both sides by 
dt / 2 , we shall obtain: 



6Tkin - 2 '^^'^^'^ 1 dt 



+ dta (3' ^]6l . (202) 



Consider the second term on the right hand side of fl202p . We know from 
ffTTTI) that 

d'j' , (203) 



so we can write the second term on the right hand side of fl202p as 



.2 - ' 1 1 



dta P {—- — ]6l = - Gd-f 51 . (204) 
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The first term on the right hand side of fl202p can be written down as 

1 1 \ / 1 1 



dt 7 57 



cos / 



-G^ dt^ 57 



l-(7')2 \I2 h) \h h 



_G7 j./i (/^')' I 

1 — (7')^ a (3' a 7' 



-G5-i' f /?'7 



7-^ (^7 + d(3 



-GS-f' 



«'(l-(7')^) 



p' 7' rfV + (1 - (7?) 



(205) 



We obtained the first equahty by noting that the definitions of the direction 
cosines yield 

while the second equahty comes from the equations of motion (11751 - I177P in 
Appendix A. 2. Finally, when we look at d(3 , we see that 

dp = d{sm J cos /) = d{^/l — (7')^ cos I) 



l^^-^Y^^Mdl = -i^-adl . (207) 
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We can now write fl205p as 
-G6^' 



«'(l-(7'n 



p' y d^' + (1 - (7?) dp') 



' /? 7 ^7 + (1 - (7 ) ) -— -a dl 



«'(l-(7')^) V ' ^ ' " ' V l-(7')^ 
= Gdlb-i 

= G6{cosJ)dl . (208) 
This enables us to write (12021) as 

6Tk^n - 2 Zk^i^ ) dt = G 5(cos J)dl-G rf(cos J) 5/ . (209) 



G 

Integration of this yields 



G j 5(cos J) rf/ = {^Tk^n - 2 IMi^^ jdt + G jsi dj' 



STk,n -2^%^^ + G COS J 61 . (210) 



Insertion of this result into (1197^ brings up the following 

5S = - ^) STkin + /i5(G cos/) + (9 + J cos J dl - 2 ^ . (211) 

Finally, getting rid of (/ cos J dl — 2(Tfcj„/G) {u/n)) by means of 
( 16T|) . we arrive at the desired formula ( 1671) . 



A. 4. Proof of formula ( 1601 ) 

Since we are now talking about unperturbed rotation, the invariable plane 
is inertial and we can employ our formulae ([S]-[7]) for the body-frame-related 
angular velocity, with angles g , J , I inserted instead of (j) , 9 , ip: 

ui = g sin J sin/ + J cos/ , (212) 
U2 = g sin J cos/ — J sin/ , (213) 
UJ3 = i + g cosJ . (214) 
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Therefrom we extract two derivatives: 

g = (ui sin / + UJ2 cos /) , 

sm J 

cos^ 6 

I cos J = UJ3 cos J (ui sin / + UJ2 cos /) , 

sm J 

summation whereof yields 

1 — cos^ J 

g + I cos J = LJ3 cos J H (ui sin I + uj2 cos I) = uj^ cos J + sin J iuji sin / + uj2 cos Z) 

sin J 

This expression, in combination with 

hoJi = G sin J sin / , I2UJ2 = G sin J cos / , I^u^ = G cos J , 
entails (|60l). 

A. 5. Andoyer's Method 

Let [xi , yi) i = 1, . . . , n , be a set of canonical variables, and let the 
Hamiltonian be some 7i(xj , yi), so that 

cbi _ dH dyj^ _ _ &H 

dt dyi ' dt dxi 

Let 2;i,2;2, • • -2:271 be some new, not necessarily canonical variables. The 
question raised by Andoyer was: How will the equations of motion look in 
these new variables? 

Making use of the Lagrange brackets. 



^1 _ / ^% _ dyi dxj \ 
^ \ du dv du dv J 
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and taking into account fl215p . we can write after some algebra 

— + [2:fc,t] + ^[2;,,z,]^ = , (216) 

a linear system wlierefrom the derivatives dzi/dt can be obtained. 
At this point, Andoyer introduces a new function given as 

i 

where K is an arbitrary function (later on it will be set zero). 
Obviously, 

The linear system fl216p may be converted into 

^Jzk /'^•^zfc (^Jzi \ dzi 

with n' = n - t. 

Let us assume now the following hypothesis: the 2n variables are split 
into two parts Pj and qj, {j = 1, .. .n), in such a way that 

Jq^ = 0, and Jp^ = Jp^iqi, t). (220) 

Hence, the previous system is split into 

#^ + ^ + 5:1^^ = 0. f^-E^^ = 0- (221) 

opk ot ^ oqj at oqk ~^ oqk dt 



i^Note that 

dH fdH dx, dH dy^ \ 

and hence, based on (j215|) . 

dH ^sp f _ dy^dxi\ _ ^ 
dzk ^ \ dzk dt dzk dt J 
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Let us now define q* = Jp. , and replace qj with these new variables. Since 
qj are functions of q* and t, we get: 

dqk ^ dql dqk 
whence the previous equations acquire the form of 
dpj _ OH' dq* _ dH' 



dt dq* ' dt dpj 



(223) 



This shows that pj are the conjugate momenta of q*. 

Let us now consider a rigid-body motion. Let 0, 6, ip be the Euler angles 
and 0, \E' be their conjugate moments. If in the definition (121 7p we set the 
function K to be zero, we shall obtain: 

J. = + e— + (224) 
oa oa oa 

where a denotes any Euler angle. By definition, a conjugate momentum is a 
partial derivative of the kinetic energy with respect to the time derivative of 
the corresponding Euler angle. Thus 

J _ dcf^dT de_dT dj^&T_ 

which we can be re-written as 

■ dT ■ dT ■ dT 

o(p od dip 

where 0a,^a,'?/'a are the virtual angular velocities with respect to a; the 
virtual motion is an instantaneous rotation about the origin, at an angular 
velocity oJq = ( ojia , uj2a , oj^a ) ■ Thus, the relation between oJq and the 
derivatives 0q,, Oa-, ipa, is the same as the relation between the angular-velocity 
vector u} and the derivatives (p, 6, ip. 

The kinetic energy is a quadratic form in the derivatives (j),9,ip. So the 
above expression shows that is the polar form of the quadratic form T 
with respect to (l)a,9a,ipa- Thus, since 

2T = hiul + hcoj + htul = Guj , (227) 



56 



we have: 



Jq, = IiUiUia + l2^2^2a + -^2^^3^^3q = G ■ UJa ■ (228) 

Having arrived to this point, Andoyer introduced the invariable plane per- 
pendicular to the angular momentum vector G, and computed the quantities 
J a in terms of the six variables I , g , h , I , G , J defined as on Fig. 3. 

Firstly, Jc = because (p , 9 , ip do not depend on G . Next, as vectors 
u}j and ujj point in the directions of the nodes on the invariable plane, 
then Jj = G ■ u}j = and Jj = G ■ u)j = . On the other hand, 
Jg = G ■ ujg = G , and Jh = G ■ uj^ = G cos / , and also Ji = G ■ ujj = 
G cos J . 

Thus, to satisfy the conditions (12201) and to get a canonical set of vari- 
ables, we must choose the new conjugate momenta as 

L = Ji = GcosJ, G = Jg, H = Jh = GcosJ. (229) 
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Figure 2: An inertial coordinate system, 81,82,83, a body-fixed frame, 
bi, b2, ba, an angular momentum-based frame and intersections of their fun- 
damental planes, denoted by i, 1, j. 
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S3 




Invariable plane 
(orthogonal to G) 



Figure 3: The same as the previous picture, but with fewer details. The 
reference coordinate system (inertial or, more generally, precessing) is con- 
stituted by axes Si , S2 , S3 . A body-fixed frame is defined by the principal 
axes bi , b2 , bs . The third frame is constituted by the angular-momentum 
vector G and a plane orthogonal thereto (the so-called invariable plane). 
The lines of nodes are denoted with i, 1, j . The attitude of the body relative 
to the reference frame is given by the Euler angles , 9 , il) . The orienta- 
tion of the invariable plane with respect to the reference frame is determined 
by the angles h and / . The inclination I is equal to the angle that the 
angular-momentum vector G makes with the reference axis S3 . The angle 
J between the invariable plane and the body equator coincides with the angle 
that G makes with the major-inerti^^xis bs of the body. The projections 
of the angular momentum toward the reference axis S3 and the body axis 
b3 are H = G cos I and L = G cos J . 
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Figure 5: The perturbed trajectory is a set of points belonging to a sequence 
of confocal instantaneous ellipses. The ellipses are not supposed to be tan- 
gent, nor even coplanar to the orbit at the intersection point. As a result, the 
physical velocity r (tangent to the trajectory at each of its points) differs 
from the Keplerian velocity g (tangent to the ellipse). 
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Figure 6: The perturbed trajectory is represented through a sequence of 
confocal instantaneous eUipses which are tangent to the trajectory at the 
intersection points, i.e., are osculating. Now, the physical velocity r (which 
is tangent to the trajectory) will coincide with the Keplerian velocity g 
(which is tangent to the ellipse) 
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